﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.DataSourcesFile;
using ESRI.ArcGIS.esriSystem;
using System.IO;
using RDotNet;
using System.Threading;
namespace ConnectWithR
{
    /// <summary>
    /// this class provides methods to apply methods and to transfer data from .NET to Statistic software R
    /// </summary>
    class RUtility
    {
        
        static int tableSeries = 0;
        public static Boolean isInitialized = false;
        public static REngine engine = null; // the connection with R
        
        /// <summary>
        /// intiate R
        /// </summary>
        /// <param name="path">the path to r.dll</param>
        public static void initialR(String path)
        {
            try
            {
                MessageForm msgForm = new MessageForm();
                msgForm.lblMessage.Text = "Loding R.......";
                new Thread((ThreadStart)delegate
                {
                    System.Windows.Forms.Application.Run(msgForm);
                }).Start();
                var envPath = Environment.GetEnvironmentVariable("PATH");
                var rBinPath = @path;
                //var rBinPath = @"C:\Program Files\R\R-2.15.1\bin\i386";
                Environment.SetEnvironmentVariable("PATH", envPath + Path.PathSeparator + rBinPath);
                engine = REngine.CreateInstance("RDotNet");
                
                engine.Initialize();
                engine.Evaluate("library(topmodel)");
                engine.Evaluate("library(raster)");
                engine.Evaluate("library(rgdal)");
                engine.Evaluate("data(huagrahuma)");
                engine.Evaluate("attach(huagrahuma)");
                engine.Evaluate("data(huagrahuma.dem)");

                RUtility.isInitialized = true;
                msgForm.Invoke((EventHandler)delegate { msgForm.Close(); });
                MessageBox.Show("Load R.dll successfully");
                TopmodelAddin.myForm.Focus();
            }
            catch (System.DllNotFoundException e)
            {
                MessageBox.Show("Can not find the R.dll in that direction");
                return;
            }
        }

        /// <summary>
        /// Get data from time series table
        /// </summary>
        /// <param name="pTable">time series tablke</param>
        /// <returns></returns>
        public static Double[] getValues(ESRI.ArcGIS.Geodatabase.ITable pTable)
        {
            Dictionary<String, Double> tsDictionary = new Dictionary<string, double>();
            List<Double> et0Aarry = new List<Double>();
            int timeIndex = pTable.FindField("TSDateTime");
            int valueIndex = pTable.FindField("TSValue");
            ESRI.ArcGIS.Geodatabase.ICursor pCursor = pTable.Search(null, true);
            ESRI.ArcGIS.Geodatabase.IRow pRow = pCursor.NextRow();
            String date = ((String)pRow.get_Value(timeIndex)).Substring(0, 10);
            Double value = 0.00;
            tsDictionary.Add(date, value);
            while (pRow != null)
            {
                date = ((String)pRow.get_Value(timeIndex)).Substring(0, 10);
                if (tsDictionary.ContainsKey(date))
                {
                    value = tsDictionary[date] + pRow.get_Value(valueIndex);
                    tsDictionary[date] = value;
                    pRow = pCursor.NextRow();
                }
                else
                {
                    value = pRow.get_Value(valueIndex);
                    tsDictionary.Add(date, value);
                    pRow = pCursor.NextRow();
                }

            }

            foreach (var element in tsDictionary)
            {
                et0Aarry.Add(element.Value);
            }

            return et0Aarry.ToArray<Double>();
        }

        /// <summary>
        /// Create simulate result table
        /// </summary>
        /// <returns>result table</returns>
        public static ITable createResultTable()
        {
            ITable pTable = null;
            IWorkspaceFactory pFact = new ShapefileWorkspaceFactory();
            IWorkspace pWorkspace = pFact.OpenFromFile("", 0);
            IFeatureWorkspace pFeatws = (IFeatureWorkspace)pWorkspace;
            String fileName = "result" + System.DateTime.Now.ToString("mmss") + ".dbf";

            if (File.Exists(fileName))
            {
                fileName = fileName.Replace(".dbf", Convert.ToString(tableSeries) + ".dbf");
            }

            UID uid = new UIDClass();
            uid.Value = "esriGeoDatabase.Object";

            IFields fields = new Fields();
            IFieldsEdit fieldsEdit = (IFieldsEdit)fields;

            IField field1 = new Field();
            IFieldEdit fieldEdit1 = (IFieldEdit)field1;

            IField field2 = new Field();
            IFieldEdit fieldEdit2 = (IFieldEdit)field2;

            IField field3 = new Field();
            IFieldEdit fieldEdit3 = (IFieldEdit)field3;

            IField field4 = new Field();
            IFieldEdit fieldEdit4 = (IFieldEdit)field4;

            IField field5 = new Field();
            IFieldEdit fieldEdit5 = (IFieldEdit)field5;

            IField field6 = new Field();
            IFieldEdit fieldEdit6 = (IFieldEdit)field6;

            IField field7 = new Field();
            IFieldEdit fieldEdit7 = (IFieldEdit)field7;

            IField field8 = new Field();
            IFieldEdit fieldEdit8 = (IFieldEdit)field8;

            IField field9 = new Field();
            IFieldEdit fieldEdit9 = (IFieldEdit)field9;

            IField field10 = new Field();
            IFieldEdit fieldEdit10 = (IFieldEdit)field10;

            IField field11 = new Field();
            IFieldEdit fieldEdit11 = (IFieldEdit)field11;

            IField field12 = new Field();
            IFieldEdit fieldEdit12 = (IFieldEdit)field12;

            IField field13 = new Field();
            IFieldEdit fieldEdit13 = (IFieldEdit)field13;

            IField field14 = new Field();
            IFieldEdit fieldEdit14 = (IFieldEdit)field14;

            IField field15 = new Field();
            IFieldEdit fieldEdit15 = (IFieldEdit)field15;


            fieldEdit1.Name_2 = "Object_id";
            fieldEdit1.Type_2 = esriFieldType.esriFieldTypeOID;
            fieldEdit1.IsNullable_2 = false;
            fieldEdit1.DefaultValue_2 = "0";
            fieldEdit1.Editable_2 = true;


            fieldEdit2.Name_2 = "FID";
            fieldEdit2.Type_2 = esriFieldType.esriFieldTypeInteger;
            fieldEdit2.IsNullable_2 = false;
            fieldEdit2.DefaultValue_2 = "0";
            fieldEdit2.Editable_2 = true;
            fieldEdit2.Length_2 = 10;

            fieldEdit3.Name_2 = "Qsim";
            fieldEdit3.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit3.IsNullable_2 = true;
            fieldEdit3.DefaultValue_2 = "0";
            fieldEdit3.Editable_2 = true;
            fieldEdit3.Length_2 = 10;

            fieldEdit4.Name_2 = "NSE";
            fieldEdit4.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit4.IsNullable_2 = true;
            fieldEdit4.DefaultValue_2 = "0";
            fieldEdit4.Editable_2 = true;
            fieldEdit4.Length_2 = 10;

            fieldEdit5.Name_2 = "qs0";
            fieldEdit5.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit5.IsNullable_2 = true;
            fieldEdit5.DefaultValue_2 = "0";
            fieldEdit5.Editable_2 = true;
            fieldEdit5.Length_2 = 10;

            fieldEdit6.Name_2 = "lnTe";
            fieldEdit6.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit6.IsNullable_2 = true;
            fieldEdit6.DefaultValue_2 = "0";
            fieldEdit6.Editable_2 = true;
            fieldEdit6.Length_2 = 10;


            fieldEdit7.Name_2 = "m";
            fieldEdit7.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit7.IsNullable_2 = true;
            fieldEdit7.DefaultValue_2 = "0";
            fieldEdit7.Editable_2 = true;
            fieldEdit7.Length_2 = 10;

            fieldEdit8.Name_2 = "Sr0";
            fieldEdit8.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit8.IsNullable_2 = true;
            fieldEdit8.DefaultValue_2 = "0";
            fieldEdit8.Editable_2 = true;
            fieldEdit8.Length_2 = 10;

            fieldEdit9.Name_2 = "Srmax";
            fieldEdit9.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit9.IsNullable_2 = true;
            fieldEdit9.DefaultValue_2 = "0";
            fieldEdit9.Editable_2 = true;
            fieldEdit9.Length_2 = 10;

            fieldEdit10.Name_2 = "td";
            fieldEdit10.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit10.IsNullable_2 = true;
            fieldEdit10.DefaultValue_2 = "0";
            fieldEdit10.Editable_2 = true;
            fieldEdit10.Length_2 = 10;

            fieldEdit11.Name_2 = "vch";
            fieldEdit11.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit11.IsNullable_2 = true;
            fieldEdit11.DefaultValue_2 = "0";
            fieldEdit11.Editable_2 = true;
            fieldEdit11.Length_2 = 10;

            fieldEdit12.Name_2 = "vr";
            fieldEdit12.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit12.IsNullable_2 = true;
            fieldEdit12.DefaultValue_2 = "0";
            fieldEdit12.Editable_2 = true;
            fieldEdit12.Length_2 = 10;

            fieldEdit13.Name_2 = "k0";
            fieldEdit13.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit13.IsNullable_2 = true;
            fieldEdit13.DefaultValue_2 = "0";
            fieldEdit13.Editable_2 = true;
            fieldEdit13.Length_2 = 10;

            fieldEdit14.Name_2 = "CD";
            fieldEdit14.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit14.IsNullable_2 = true;
            fieldEdit14.DefaultValue_2 = "0";
            fieldEdit14.Editable_2 = true;
            fieldEdit14.Length_2 = 10;

            fieldEdit15.Name_2 = "dt";
            fieldEdit15.Type_2 = esriFieldType.esriFieldTypeDouble;
            fieldEdit15.IsNullable_2 = true;
            fieldEdit15.DefaultValue_2 = "0";
            fieldEdit15.Editable_2 = true;
            fieldEdit15.Length_2 = 10;

            fieldsEdit.FieldCount_2 = 15;

            fieldsEdit.set_Field(0, fieldEdit1);
            fieldsEdit.set_Field(1, fieldEdit2);
            fieldsEdit.set_Field(2, fieldEdit3);
            fieldsEdit.set_Field(3, fieldEdit4);
            fieldsEdit.set_Field(4, fieldEdit5);
            fieldsEdit.set_Field(5, fieldEdit6);
            fieldsEdit.set_Field(6, fieldEdit7);
            fieldsEdit.set_Field(7, fieldEdit8);
            fieldsEdit.set_Field(8, fieldEdit9);
            fieldsEdit.set_Field(9, fieldEdit10);
            fieldsEdit.set_Field(10, fieldEdit11);
            fieldsEdit.set_Field(11, fieldEdit12);
            fieldsEdit.set_Field(12, fieldEdit13);
            fieldsEdit.set_Field(13, fieldEdit14);
            fieldsEdit.set_Field(14, fieldEdit15);


            fields = (IFields)fieldsEdit;

            IFieldChecker fieldChecker = new FieldChecker();
            IEnumFieldError enumFieldError = null;
            IFields validateFields = null;
            fieldChecker.ValidateWorkspace = pWorkspace;
            fieldChecker.Validate(fields, out enumFieldError, out validateFields);

            pTable = pFeatws.CreateTable(fileName, validateFields, uid, null, "");

            MessageBox.Show("result table created");
            return pTable;
            
        }

        private static ITable insertRow(ITable pTable, int fid, Double Qsim, Double NSE, Double[] parameter,ICursor insertCursor, IRowBuffer iRowBuf)
        {

            int oidIndex = pTable.FindField(pTable.OIDFieldName);
            int fidIndex = pTable.FindField("FID_");
            int QsimIndex = pTable.FindField("Qsim");
            int NSEIndex = pTable.FindField("NSE");
            int qs0Index = pTable.FindField("qs0");
            int lnTeIndex = pTable.FindField("lnTe");
            int mIndex = pTable.FindField("m");
            int Sr0Index = pTable.FindField("Sr0");
            int SrmaxIndex = pTable.FindField("Srmax");
            int tdIndex = pTable.FindField("td");
            int vchIndex = pTable.FindField("vch");
            int vrIndex = pTable.FindField("vr");
            int k0Index = pTable.FindField("k0");
            int CDIndex = pTable.FindField("CD");
            int dtIndex = pTable.FindField("dt");

            iRowBuf.set_Value(fidIndex, fid);
            iRowBuf.set_Value(QsimIndex, Qsim);
            iRowBuf.set_Value(NSEIndex, NSE);
            iRowBuf.set_Value(qs0Index, parameter[0]);
            iRowBuf.set_Value(lnTeIndex, parameter[1]);
            iRowBuf.set_Value(mIndex, parameter[2]);
            iRowBuf.set_Value(Sr0Index, parameter[3]);
            iRowBuf.set_Value(SrmaxIndex, parameter[4]);
            iRowBuf.set_Value(tdIndex, parameter[5]);
            iRowBuf.set_Value(vchIndex, parameter[6]);
            iRowBuf.set_Value(vrIndex, parameter[7]);
            iRowBuf.set_Value(k0Index, parameter[8]);
            iRowBuf.set_Value(CDIndex, parameter[9]);
            iRowBuf.set_Value(dtIndex, parameter[10]);

            insertCursor.InsertRow(iRowBuf);
            insertCursor.Flush();
            return pTable;

        }

        /// <summary>
        /// calculate simulated discharge
        /// </summary>
        /// <param name="parameters">parameter array</param>
        /// <param name="DEMPath">path of DEM</param>
        /// <param name="PrecipitationValues">precipitation array</param>
        /// <param name="ET0Values">evapotranspiration array</param>
        /// <param name="drainageArea">acreage of search area</param>
        /// <returns></returns>
        public static ITable calculateCustomDischarge(Double[] parameters, String DEMPath, Double[] PrecipitationValues, Double[] ET0Values, Double drainageArea)
        {
            if (PrecipitationValues.Length != ET0Values.Length)
            {
                MessageBox.Show("Precipitation data and evapotranspiration data must have the same length!");
                return null;
            }
            else if (RUtility.engine == null)
            {
                MessageBox.Show("please select r.dll directory first");
                return null;
            }
            else
            {

                MessageForm msgForm = new MessageForm();
                msgForm.lblMessage.Text = "Calculating...";
                new Thread((ThreadStart)delegate
                {
                    System.Windows.Forms.Application.Run(msgForm);
                }).Start();
                NumericVector ET0 = engine.CreateNumericVector(ET0Values);
                engine.SetSymbol("ET0Cuntom", ET0);
                NumericVector rain = engine.CreateNumericVector(PrecipitationValues);
                engine.SetSymbol("rainCuntom", rain);
                NumericVector para = engine.CreateNumericVector(parameters);
                engine.SetSymbol("parametersCuntom", para);
                DEMPath = DEMPath.Replace("\\", "\\\\");
                String sentence = "demCuntom <- raster(\"" + DEMPath + "\")";
                engine.Evaluate(sentence);
                sentence = "demCuntom <- as.matrix(demCuntom)";
                engine.Evaluate(sentence);
                sentence = "result1 <- demCuntom[rowSums(is.na(demCuntom))!= dim(demCuntom)[2],]";
                engine.Evaluate(sentence);
                sentence = "result2 <- result1[,colSums(is.na(result1))!= dim(result1)[1]]";
                engine.Evaluate(sentence);
                sentence = "topix <- topidx(result2, resolution=10)";
                engine.Evaluate(sentence);
                sentence = "topix_after <- make.classes(topix$atb,20)";
                engine.Evaluate(sentence);
                sentence = "Qsim <- topmodel(parametersCuntom, topix_after, delay, rainCuntom, ET0Cuntom)";
                engine.Evaluate(sentence);
                Double[] qsim = engine.GetSymbol("Qsim").AsNumeric().ToArray<Double>();

                //Transform discharge unit from m/day to cfs
                for (int i = 0; i < qsim.Length; i++)
                {
                    qsim[i] = qsim[i] / 86400 * drainageArea * 35.3147;
                }
                msgForm.Invoke((EventHandler)delegate { msgForm.Close(); });


                ITable resultTable = createResultTable();
                ICursor insertCursor = resultTable.Insert(true);
                IRowBuffer iRowBuf = resultTable.CreateRowBuffer();
                for (int i = 0; i < qsim.Length; i++)
                {
                    insertRow(resultTable, i, qsim[i], -1, parameters, insertCursor, iRowBuf);
                }
                return resultTable;
            }
        }

        /// <summary>
        /// calcualte Nash Sutcliffe efficiencies
        /// </summary>
        /// <param name="parameters">parameter array</param>
        /// <param name="DEMPath">path of DEM</param>
        /// <param name="PrecipitationValues">precipitation array</param>
        /// <param name="ET0Values">evapotranspiration array</param>
        /// <param name="dischargeValues">discharge array</param>
        /// <param name="drainageArea">acreage of research area</param>
        /// <returns></returns>
        public static ITable calculateCustomNSE(Double[] parameters, String DEMPath, Double[] PrecipitationValues, Double[] ET0Values, Double[] dischargeValues, Double drainageArea)
        {
            if (PrecipitationValues.Length != ET0Values.Length || PrecipitationValues.Length!=dischargeValues.Length || dischargeValues.Length != ET0Values.Length)
            {
                MessageBox.Show("Precipitation data and evapotranspiration data must have the same length!");
                return null;
            }
            else if (RUtility.engine == null)
            {
                MessageBox.Show("please select r.dll directory first");
                return null;
            }
            else
            {
                MessageForm msgForm = new MessageForm();
                msgForm.lblMessage.Text = "Calculating...";
                new Thread((ThreadStart)delegate
                {
                    System.Windows.Forms.Application.Run(msgForm);
                }).Start();

                NumericVector ET0 = engine.CreateNumericVector(ET0Values);
                engine.SetSymbol("ET0Custom_NSE", ET0);
                NumericVector rain = engine.CreateNumericVector(PrecipitationValues);
                engine.SetSymbol("rainCustom_NSE", rain);
                NumericVector discharge = engine.CreateNumericVector(dischargeValues);
                engine.SetSymbol("dischargeCustom_NSE", discharge);
                NumericVector para = engine.CreateNumericVector(parameters);
                engine.SetSymbol("parametersCustom_NSE", para);

                DEMPath = DEMPath.Replace("\\", "\\\\");
                String sentence = "demCuntom <- raster(\"" + DEMPath + "\")";
                engine.Evaluate(sentence);
                sentence = "demCuntom <- as.matrix(demCuntom)";
                engine.Evaluate(sentence);
                sentence = "result1 <- demCuntom[rowSums(is.na(demCuntom))!= dim(demCuntom)[2],]";
                engine.Evaluate(sentence);
                sentence = "result2 <- result1[,colSums(is.na(result1))!= dim(result1)[1]]";
                engine.Evaluate(sentence);
                sentence = "topix <- topidx(result2, resolution=10)";
                engine.Evaluate(sentence);
                sentence = "topix_after__NSE <- make.classes(topix$atb,20)";
                engine.Evaluate(sentence);
                sentence = "Qsim_NSE <- topmodel(parametersCustom_NSE, topix_after_NSE, delay, rainCustom_NSE, ET0Custom_NSE )";
                engine.Evaluate(sentence);
                sentence = "nse_NSE <- topmodel(parametersCustom_NSE, topix_after_NSE, delay, rainCustom_NSE, ET0Custom_NSE, Qobs = dischargeCustom_NSE)";
                engine.Evaluate(sentence);
                Double[] qsim = engine.GetSymbol("Qsim_NSE").AsNumeric().ToArray<Double>();
                Double[] nse = engine.GetSymbol("nse_NSE").AsNumeric().ToArray<Double>();

                //Transform discharge unit from m/day to cfs
                for (int i = 0; i < qsim.Length; i++)
                {
                    qsim[i] = qsim[i] / 86400 * drainageArea * 35.3147;
                }
                msgForm.Invoke((EventHandler)delegate { msgForm.Close(); });


                ITable resultTable = createResultTable();
                ICursor insertCursor = resultTable.Insert(true);
                IRowBuffer iRowBuf = resultTable.CreateRowBuffer();
                for (int i = 0; i < qsim.Length; i++)
                {
                    insertRow(resultTable, i, qsim[i], nse[0], parameters, insertCursor, iRowBuf);
                }

                return resultTable;
            }
        }

        /// <summary>
        /// calcualte Nash Sutcliffe efficiencies array by a Monte carlo sampling
        /// </summary>
        /// <param name="parameters">parameter array</param>
        /// <param name="DEMPath">path of DEM</param>
        /// <param name="PrecipitationValues">precipitation array</param>
        /// <param name="ET0Values">evapotranspiration array</param>
        /// <param name="dischargeValues">discharge array</param>
        /// <param name="drainageArea">acreage of research area</param>
        /// <returns></returns>
        public static ITable calculateMonteCarlo(Double[] parameters, String DEMPath, Double[] PrecipitationValues, Double[] ET0Values, Double[] dischargeValues)
        {
            if (PrecipitationValues.Length != ET0Values.Length || PrecipitationValues.Length != dischargeValues.Length || dischargeValues.Length != ET0Values.Length)
            {
                MessageBox.Show("Precipitation data and evapotranspiration data must have the same length!");
                return null;
            }
            else if (RUtility.engine == null)
            {
                MessageBox.Show("please select r.dll directory first");
                return null;
            }
            else
            {
                MessageForm msgForm = new MessageForm();
                msgForm.lblMessage.Text = "Calculating...";
                new Thread((ThreadStart)delegate
                {
                    System.Windows.Forms.Application.Run(msgForm);
                }).Start();

                Double[] qs0Array = { parameters[0] };
                Double[] inTeArray = { parameters[1] };
                Double[] mArray = { parameters[2] };
                Double[] Sr0Array = { parameters[3] };
                Double[] SrmaxArray = { parameters[4] };
                Double[] tdArray = { parameters[5] };
                Double[] vrArray = { parameters[6] };
                Double[] k0Array = { parameters[7] };
                Double[] cdArray = { parameters[8] };
                Double[] runsArray = { parameters[9] };

                String sentence = "rm(list=ls())";
                engine.Evaluate(sentence);
                NumericVector qs0ArrayVevtor = engine.CreateNumericVector(qs0Array);
                engine.SetSymbol("qs0max", qs0ArrayVevtor);
                NumericVector inTeArrayVevtor = engine.CreateNumericVector(inTeArray);
                engine.SetSymbol("inTemax", inTeArrayVevtor);
                NumericVector mArrayVevtor = engine.CreateNumericVector(mArray);
                engine.SetSymbol("mmax", mArrayVevtor);
                NumericVector Sr0ArrayVevtor = engine.CreateNumericVector(Sr0Array);
                engine.SetSymbol("Sr0max", Sr0ArrayVevtor);
                NumericVector SrmaxArrayVevtor = engine.CreateNumericVector(SrmaxArray);
                engine.SetSymbol("Srmaxmax", SrmaxArrayVevtor);
                NumericVector tdArrayVevtor = engine.CreateNumericVector(tdArray);
                engine.SetSymbol("tdmax", tdArrayVevtor);
                NumericVector vrArrayVevtor = engine.CreateNumericVector(vrArray);
                engine.SetSymbol("vrmax", vrArrayVevtor);
                NumericVector k0ArrayVevtor = engine.CreateNumericVector(k0Array);
                engine.SetSymbol("k0max", k0ArrayVevtor);
                NumericVector cdArrayVevtor = engine.CreateNumericVector(cdArray);
                engine.SetSymbol("cdmax", cdArrayVevtor);
                NumericVector runsArrayVevtor = engine.CreateNumericVector(runsArray);
                engine.SetSymbol("runs", runsArrayVevtor);


                NumericVector ET0 = engine.CreateNumericVector(ET0Values);
                engine.SetSymbol("ET0Custom_Monte", ET0);
                NumericVector rain = engine.CreateNumericVector(PrecipitationValues);
                engine.SetSymbol("rainCustom_Monte", rain);

                NumericVector discharge = engine.CreateNumericVector(dischargeValues);
                engine.SetSymbol("dischargeCustom_Monte", discharge);

                DEMPath = DEMPath.Replace("\\", "\\\\");
                sentence = "demCuntom <- raster(\"" + DEMPath + "\")";
                engine.Evaluate(sentence);
                sentence = "demCuntom <- as.matrix(demCuntom)";
                engine.Evaluate(sentence);
                sentence = "result1 <- demCuntom[rowSums(is.na(demCuntom))!= dim(demCuntom)[2],]";
                engine.Evaluate(sentence);
                sentence = "result2 <- result1[,colSums(is.na(result1))!= dim(result1)[1]]";
                engine.Evaluate(sentence);
                sentence = "topix <- topidx(result2, resolution=10)";
                engine.Evaluate(sentence);
                sentence = "topix_after_Monte <- make.classes(topix$atb,20)";
                engine.Evaluate(sentence);
                sentence = "qs0 <- runif(runs, min = 0, max = qs0max)";
                engine.Evaluate(sentence);
                sentence = "inTe <- runif(runs, min = -2, max = inTemax)";
                engine.Evaluate(sentence);
                sentence = "m <- runif(runs, min = 0, max = mmax)";
                engine.Evaluate(sentence);
                sentence = "Sr0 <- runif(runs, min = 0, max = Sr0max)";
                engine.Evaluate(sentence);
                sentence = "Srmax <- runif(runs, min = 0, max = Srmaxmax)";
                engine.Evaluate(sentence);
                sentence = "td <- runif(runs, min = 0, max = tdmax)";
                engine.Evaluate(sentence);
                sentence = "vch <- 1000";
                engine.Evaluate(sentence);
                sentence = "vr <- runif(runs, min = 0, max = vrmax)";
                engine.Evaluate(sentence);
                sentence = "k0 <- runif(runs, min = 0, max = k0max)";
                engine.Evaluate(sentence);
                sentence = "cd <- runif(runs, min = 0, max = cdmax)";
                engine.Evaluate(sentence);
                sentence = "dt <- 1";
                engine.Evaluate(sentence);
                sentence = "parametersCustom_Monte <- cbind(qs0, inTe, m, Sr0, Srmax, td, vch, vr, k0, cd, dt)";
                engine.Evaluate(sentence);
                String sentence1 = "nse_Monte <- topmodel(parametersCustom_Monte, topix_after_Monte, delay, rainCustom_Monte, ET0Custom_Monte, Qobs = dischargeCustom_Monte)";
                engine.Evaluate(sentence1);

                NumericMatrix parametersCustom = engine.GetSymbol("parametersCustom_Monte").AsNumericMatrix();
                Double[,] parametersCollection = new Double[parametersCustom.RowCount, parametersCustom.ColumnCount];
                parametersCustom.CopyTo(parametersCollection, parametersCustom.RowCount, parametersCustom.ColumnCount);
                Double[] nse = engine.GetSymbol("nse_Monte").AsNumeric().ToArray<Double>();

                msgForm.Invoke((EventHandler)delegate { msgForm.Close(); });


                ITable resultTable = createResultTable();
                ICursor insertCursor = resultTable.Insert(true);
                IRowBuffer iRowBuf = resultTable.CreateRowBuffer();
                for (int i = 0; i < nse.Length; i++)
                {
                    List<Double> parameterList  = new List<double>();
                    for(int j = 0 ; j<parametersCustom.ColumnCount; j++)
                    {
                        parameterList.Add(parametersCollection[i,j]);
                    }
                    insertRow(resultTable, i, -1, nse[i], parameterList.ToArray<Double>(), insertCursor, iRowBuf);
                }

                sentence = null;
                insertCursor = null;
                iRowBuf = null;
                return resultTable;
            }
        }
    }
}
